Models Cholesterol_prior_sens and asprin_logistic_random_1_priors_sensitivity originally failed with Attempt to redefine node errors. I got them to run by moving the troublesome lines to data blocks. These chunks are marked with eval=TRUE.

TO DO:

library('rjags')
library('coda')
library('jagshelper')
library('jagsUI')
library('mcmcplots')  # caterplot
library('dplyr')

Aspirin

asprin_normal_fixed_1.txt

This is the model from Figure 4.1 (with a few extra digits of precision in the data).

model_code <- "model 
{ 
  for (i in 1:Nstud)
    {
      P[i] <- 1/V[i]
      y[i] ~ dnorm(d, P[i])
    }
    d ~ dnorm(0, 1.0E-5)
    OR <- exp(d)
}" %>% textConnection

# DATA
data <- list(
  y=c(.3289011, .3845458, .2195622, .2222206, .2254672, -.1246363, .1109658), 
  V=c(.0388957,.0411673,.0204915,.0647646,.0351996,.0096167,.0015062), 
  Nstud=7
)

# INITIAL VALUES: This is a list of lists, one per chain.
initial_values <- list(
  list(d=0)
)

results <- jags(
  data = data,
  inits = initial_values,
  parameters.to.save = c("OR", "d"),
  model.file = model_code,
  n.chains = length(initial_values),
  n.adapt = 100,
  n.iter = 50000,
  n.burnin = 20000,
  n.thin = 2
)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 7
##    Unobserved stochastic nodes: 1
##    Total graph size: 27
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 1 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 20000 iterations x 1 chains 
##  
## 
## Sampling from joint posterior, 30000 iterations x 1 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.
summary(results)
##                mean         sd        2.5%         25%        50%        75%
## OR        1.1154942 0.03697400  1.04527705  1.09023746  1.1149566  1.1404389
## d         0.1087484 0.03314299  0.04428197  0.08639553  0.1088155  0.1314132
## deviance -3.7395092 1.41023676 -4.73853030 -4.63900994 -4.2773680 -3.4146738
##             97.5% Rhat n.eff overlap0      f
## OR       1.189913   NA     1        0 1.0000
## d        0.173880   NA     1        0 0.9996
## deviance 0.288420   NA     1        1 0.9704
plot(results)

asprin_normal_random_1.txt

This is the model from Figure 4.3, except it uses lower case for the y[] variable vector.

model_code <- "model 
{ 
  for (i in 1:Nstud)
    {
      P[i] <- 1/V[i]
      y[i] ~ dnorm(delta[i], P[i])
      delta[i] ~ dnorm(d, prec) 
    }   
    d ~ dnorm(0, 1.0E-5)
    OR <- exp(d)
    tau ~ dunif(0,10)
    tau.sq <- tau*tau
    prec <- 1/(tau.sq)
}" %>% textConnection

# DATA  (the same as for the fixed effects model above)
data <- list(
  y=c(.3289011, .3845458, .2195622, .2222206, .2254672, -.1246363, .1109658), 
  V=c(.0388957,.0411673,.0204915,.0647646,.0351996,.0096167,.0015062), 
  Nstud=7
)

# INITIAL VALUES: This is a list of lists, one per chain.
initial_values <- list(
  list(d = 0, tau = 1, delta = c(0,0,0,0,0,0,0))
)

results <- jags(
  data = data,
  inits = initial_values,
  parameters.to.save = c("OR", "d", "prec", "tau.sq"),
  model.file = model_code,
  n.chains = length(initial_values),
  n.adapt = 100,
  n.iter = 50000,
  n.burnin = 20000,
  n.thin = 2
)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 7
##    Unobserved stochastic nodes: 9
##    Total graph size: 38
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 1 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 20000 iterations x 1 chains 
##  
## 
## Sampling from joint posterior, 30000 iterations x 1 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.
summary(results)
##                   mean           sd          2.5%           25%         50%
## OR        1.160938e+00 1.122747e-01  9.791676e-01   1.091816793  1.14621597
## d         1.447896e-01 9.327927e-02 -2.105249e-02   0.087843091  0.13646606
## prec      1.199015e+04 5.899705e+05  5.368759e+00  22.382149023 51.48911368
## tau.sq    3.712591e-02 5.701788e-02  7.957887e-05   0.006451449  0.01942158
## deviance -7.852934e+00 3.351394e+00 -1.306269e+01 -10.498470035 -8.27432876
##                   75%         97.5% Rhat n.eff overlap0         f
## OR         1.21432738     1.4237089   NA     1        0 1.0000000
## d          0.19419032     0.3532654   NA     1        1 0.9621333
## prec     155.00394007 12566.1511858   NA     1        0 1.0000000
## tau.sq     0.04467846     0.1862628   NA     1        0 1.0000000
## deviance  -5.48017083    -0.7749141   NA     1        0 0.9838000
plot(results)

asprin_logistic_random_1.txt

This is the model from Figure 4.5 (except that the line OR <- exp(d) is in a different place).

model_code <- "model
{
  for( i in 1 : Nstud ) {
    rA[i] ~ dbin(pA[i], nA[i])
    rB[i] ~ dbin(pB[i], nB[i])
    logit(pA[i]) <- mu[i]
    logit(pB[i]) <- mu[i] + delta[i]
    mu[i] ~ dnorm(0.0,1.0E-5)
    delta[i] ~ dnorm(d, prec)
  }
  OR <- exp(d)   # !!! this line comes last in Figure 4.5
  d ~ dnorm(0.0,1.0E-6)
  tau ~ dunif(0,10)
  tau.sq <- tau*tau
  prec <- 1/(tau.sq)
}" %>% textConnection

# DATA
data <- list(
  rA=c(49,44,102,32,85,246,1570), 
  rB=c(67,64,126,38,52,219,1720), 
  nA=c(615,758,832,317,810,2267,8587), 
  nB=c(624,771,850,309,406,2257,8600), 
  Nstud=7
)


# INITIAL VALUES: This is a list of lists, one per chain.
initial_values <- list(
  list(d=0, tau=1, delta=c(0,0,0,0,0,0,0), mu=c(0,0,0,0,0,0,0))
)

results <- jags(
  data = data,
  inits = initial_values,
  parameters.to.save = c("OR", "d", "prec", "tau.sq"),
  model.file = model_code,
  n.chains = length(initial_values),
  n.adapt = 100,
  n.iter = 50000,
  n.burnin = 20000,
  n.thin = 2
)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 14
##    Unobserved stochastic nodes: 16
##    Total graph size: 74
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 1 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 20000 iterations x 1 chains 
##  
## 
## Sampling from joint posterior, 30000 iterations x 1 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.
summary(results)
##                  mean           sd          2.5%          25%          50%
## OR       1.166672e+00 1.187119e-01  0.9813849004 1.095177e+00   1.15035462
## d        1.493962e-01 9.621368e-02 -0.0187905424 9.091595e-02   0.14007026
## prec     1.492055e+04 5.250138e+05  5.4773630417 2.195923e+01  48.64527803
## tau.sq   3.808501e-02 6.082000e-02  0.0001689156 7.508704e-03   0.02055698
## deviance 1.044143e+02 5.056343e+00 96.3293298236 1.007238e+02 103.84752807
##                   75%        97.5% Rhat n.eff overlap0         f
## OR         1.22119445    1.4374136   NA     1        0 1.0000000
## d          0.19982943    0.3628454   NA     1        1 0.9635333
## prec     133.17880903 5920.1166164   NA     1        0 1.0000000
## tau.sq     0.04553893    0.1825696   NA     1        0 1.0000000
## deviance 107.44671303  115.9898590   NA     1        0 1.0000000
plot(results)

asprin_logistic_random1_predict.txt

model_code <- "model
{
  for( i in 1 : Nstud ) {
    rA[i] ~ dbin(pA[i], nA[i])
    rB[i] ~ dbin(pB[i], nB[i])
    logit(pA[i]) <- mu[i]
    logit(pB[i]) <- mu[i] + delta[i]
    mu[i] ~ dnorm(0.0,1.0E-5)
    delta[i] ~ dnorm(d, prec)
  }
  
  OR <- exp(d)
  d ~ dnorm(0.0,1.0E-6)
  tau~dunif(0,10)
  tau.sq<-tau*tau
  prec<-1/(tau.sq)
  
  d.new ~ dnorm(d, prec)
  OR.new <-exp(d.new)
}" %>% textConnection

# DATA
data <- list(
  rA=c(49,44,102,32,85,246,1570), 
  rB=c(67,64,126,38,52,219,1720), 
  nA=c(615,758,832,317,810,2267,8587), 
  nB=c(624,771,850,309,406,2257,8600), 
  Nstud=7
)


# INITIAL VALUES: This is a list of lists, one per chain.
initial_values <- list(
  list(d=0, tau =1, delta=c(0,0,0,0,0,0,0), mu=c(0,0,0,0,0,0,0))
)

results <- jags(
  data = data,
  inits = initial_values,
  parameters.to.save = c("OR", "OR.new", "d.new", "delta", "tau.sq"),
  model.file = model_code,
  n.chains = length(initial_values),
  n.adapt = 100,
  n.iter = 50000,
  n.burnin = 20000,
  n.thin = 2
)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 14
##    Unobserved stochastic nodes: 17
##    Total graph size: 76
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 1 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 20000 iterations x 1 chains 
##  
## 
## Sampling from joint posterior, 30000 iterations x 1 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.
summary(results)
##                  mean         sd         2.5%           25%          50%
## OR         1.16482673 0.11727540  0.978232230   1.093407115   1.15076623
## OR.new     1.18844654 0.30154239  0.737259168   1.040355891   1.14289962
## d.new      0.14666203 0.22323556 -0.304815798   0.039562857   0.13356856
## delta[1]   0.20910410 0.14090445 -0.031117108   0.111941872   0.19040952
## delta[2]   0.22763201 0.14769461 -0.015417530   0.122630883   0.20781802
## delta[3]   0.17825449 0.11126175 -0.019969587   0.102755952   0.16876713
## delta[4]   0.16310130 0.14905091 -0.119012662   0.072870934   0.15053163
## delta[5]   0.17095064 0.12942249 -0.070930538   0.087314807   0.15946323
## delta[6]  -0.02882947 0.10186744 -0.238741923  -0.098195039  -0.02536030
## delta[7]   0.11311926 0.03778339  0.037838460   0.087995759   0.11337927
## tau.sq     0.04009716 0.06411076  0.000241616   0.008282165   0.02212963
## deviance 104.33057397 4.97948307 96.211258608 100.750212169 103.83358607
##                   75%       97.5% Rhat n.eff overlap0         f
## OR         1.21896888   1.4333353   NA     1        0 1.0000000
## OR.new     1.27782385   1.8985090   NA     1        0 1.0000000
## d.new      0.24515851   0.6410688   NA     1        1 0.8200000
## delta[1]   0.29308486   0.5260593   NA     1        1 0.9574000
## delta[2]   0.31808150   0.5573990   NA     1        1 0.9666667
## delta[3]   0.24675171   0.4168797   NA     1        1 0.9593333
## delta[4]   0.24749593   0.4904861   NA     1        1 0.8915333
## delta[5]   0.24855173   0.4524951   NA     1        1 0.9276667
## delta[6]   0.04576045   0.1522363   NA     1        1 0.5945333
## delta[7]   0.13863053   0.1865140   NA     1        0 0.9986667
## tau.sq     0.04677193   0.1944284   NA     1        0 1.0000000
## deviance 107.32915712 115.4040223   NA     1        0 1.0000000
plot(results)

asprin_logistic_random_1_priors_sensitivity

RUNTIME ERROR:
Compilation error on line 13.
Attempt to redefine node rBx[1,1]

I moved the offending lines into a data block, as recommended on stackoverflow.

Original WinBUGS code:

model {
    
  ## replicating the data 2 times to fit to 3 models in total where j indexes alternative prior distributions
  
  for (j in 1:3)
    {
    
      for( i in 1 : Nstud )
      {
        rAx[j,i] ~ dbin(pA[j,i], nA[i])
        rAx[j,i] <- rA[i]
        rBx[j,i] ~ dbin(pB[j,i], nB[i])
        rBx[j,i] <- rB[i]
        logit(pA[j,i]) <- mu[j,i]
        logit(pB[j,i]) <- mu[j,i] + delta[j,i]
        mu[j,i] ~  dnorm(0.0,1.0E-5)
        delta[j,i] ~ dnorm(d[j], prec[j])
      }
    }
  
  for (j in 1:3)
    {
      OR[j] <- exp(d[j])
      d[j] ~ dnorm(0.0,1.0E-6)
      m[j]~dnorm(0.0,1.0E-6)
    }

    
  # Setting three different priors for the between study variance
  
  prec[1] <- 1/tau.sq[1]
  tau.sq[1] <- tau[1]*tau[1]
  tau[1] ~ dnorm(0,1.0E-6)I(0,)
  
  prec[2] ~ dgamma(0.001,0.001)
  tau.sq[2] <- 1/prec[2]
  tau[2] <- sqrt(tau.sq[2])
  
  tau[3] ~ dunif(0,10)
  tau.sq[3] <- tau[3]*tau[3]
  prec[3] <- 1/tau.sq[3]
    
}
model_code <- "

data {
  for (j in 1:3) {
    for( i in 1 : Nstud ) {
      rAx[j,i] <- rA[i]
      rBx[j,i] <- rB[i]
    }
  }
}

model{
    
  ## replicating the data 2 times to fit to 3 models in total where j indexes alternative prior distributions
  
  for (j in 1:3) {
    for( i in 1 : Nstud ) {
      rAx[j,i] ~ dbin(pA[j,i], nA[i])
      # rAx[j,i] <- rA[i]
      rBx[j,i] ~ dbin(pB[j,i], nB[i])
      # rBx[j,i] <- rB[i]
      logit(pA[j,i]) <- mu[j,i]
      logit(pB[j,i]) <- mu[j,i] + delta[j,i]
      mu[j,i] ~  dnorm(0.0,1.0E-5)
      delta[j,i] ~ dnorm(d[j], prec[j])
    }
  }
  
  for (j in 1:3){
    OR[j] <- exp(d[j])
    d[j] ~ dnorm(0.0,1.0E-6)
    m[j]~dnorm(0.0,1.0E-6)
  }

    
  # Setting three different priors for the between study variance
  
  prec[1] <- 1/tau.sq[1]
  tau.sq[1] <- tau[1]*tau[1]
  tau[1] ~ dnorm(0,1.0E-6)I(0,)
  
  prec[2] ~ dgamma(0.001,0.001)
  tau.sq[2] <- 1/prec[2]
  tau[2] <- sqrt(tau.sq[2])
  
  tau[3] ~ dunif(0,10)
  tau.sq[3] <- tau[3]*tau[3]
  prec[3] <- 1/tau.sq[3]
    
}
" %>% textConnection

# DATA
data <- list(
  rA=c(49,44,102,32,85,246,1570), 
  rB=c(67,64,126,38,52,219,1720), 
  nA=c(615,758,832,317,810,2267,8587), 
  nB=c(624,771,850,309,406,2257,8600), 
  Nstud=7
)


# INITIAL VALUES: This is a list of lists, one per chain.
initial_values <- list(
  # WinBUGS code has extra closing parentheses (???)
  list(
    prec=c(NA,0.1,NA),  
    tau=c(0.01, NA, 0.01), 
    mu = structure(
      .Data=c(0,0,0,0,0,0,0,  0,0,0,0,0,0,0,  0,0,0,0,0,0,0),
      .Dim=c(3,7)),
    delta = structure(
      .Data=c(0,0,0,0,0,0,0,   0,0,0,0,0,0,0,   0,0,0,0,0,0,0),
      .Dim=c(3,7)), 
    d =c(0.1,0.1,0.1)
  )

)

results <- jags(
  data = data,
  inits = initial_values,
  parameters.to.save = c("OR", "d"),
  model.file = model_code,
  n.chains = length(initial_values),
  n.adapt = 100,
  n.iter = 50000,
  n.burnin = 20000
)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling data graph
##    Resolving undeclared variables
##    Allocating nodes
##    Initializing
##    Reading data back into data table
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 42
##    Unobserved stochastic nodes: 51
##    Total graph size: 200
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 1 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 20000 iterations x 1 chains 
##  
## 
## Sampling from joint posterior, 30000 iterations x 1 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.
summary(results)
##                 mean         sd          2.5%          25%        50%
## OR[1]      1.1662222 0.11739451   0.981212714   1.09437789   1.149769
## OR[2]      1.1526115 0.09525176   0.997237644   1.09143046   1.141165
## OR[3]      1.1650587 0.11863985   0.976503633   1.09352852   1.149952
## d[1]       0.1489909 0.09671732  -0.018966010   0.09018606   0.139561
## d[2]       0.1387560 0.08024680  -0.002766179   0.08748919   0.132050
## d[3]       0.1479067 0.09766691  -0.023776808   0.08940964   0.139720
## deviance 313.4448436 8.64787962 297.834584100 307.38539084 312.924119
##                  75%       97.5% Rhat n.eff overlap0         f
## OR[1]      1.2216751   1.4400013   NA     1        0 1.0000000
## OR[2]      1.2004897   1.3714749   NA     1        0 1.0000000
## OR[3]      1.2194936   1.4332097   NA     1        0 1.0000000
## d[1]       0.2002229   0.3646440   NA     1        1 0.9625333
## d[2]       0.1827296   0.3158868   NA     1        1 0.9730000
## d[3]       0.1984357   0.3599164   NA     1        1 0.9605000
## deviance 319.0349227 331.7744037   NA     1        0 1.0000000
plot(results)

Cholesterol

chol_model.txt, chol_data.txt

Error parsing model file:
syntax error on line 5 near "var"

I’m not sure why this was in the Supplemental Materials. It seems to be exactly the same as Cholesterol_fixed.txt below (Model 4.1 in the Exercises), except it uses the name ‘var’ instead of ‘V’. It turns out that ‘var’ is problematic in JAGS; maybe it is a reserved word or something? I changed it to V (as in Cholesterol_fixed.txt), which seems to work.

model_code <- "model 
{ 
  for (i in 1:k)
  {
    precision[i] <- 1/V[i]
    logor[i] ~ dnorm(theta, precision[i])
  }
    theta ~ dnorm(0, 1.0E-5)
    OR <- exp(theta)
}
" %>% textConnection

# DATA
data <- list(k=34)

data_df <- "logor   V
0.4893853   0.1987757
-0.0683811  0.054427
-0.4876372  0.0730886
-0.1630902  0.1108992
1.469676    1.186739
0.0057637   0.0608936
-0.1496961  0.0799553
-1.94591    2.351288
-0.0470209  0.0194635
-0.3824268  0.0552552
-0.5508629  0.0615112
-0.0984401  0.0510244
0.2939126   0.1708592
-0.4883528  1.582492
-0.4354064  0.0030637
-0.2885568  0.0496847
0.4696847   0.0716345
0.2692525   0.0101539
-1.563976   0.8819983
-0.3517397  0.3733933
-0.0480808  0.0298812
-1.109251   2.687944
0.0129744   0.0469846
-0.4124823  0.0383304
-0.0233855  0.0200697
0.0803459   0.0082181
-0.2820902  0.0420423
-1.098612   2.707904
0.521464    2.696409
-0.2501722  0.477032
1.025401    0.3643915
-0.7528253  0.0676265
-1.580711   2.272708
0.5030903   0.1426682" %>% read.delim(text=., sep="\t")

for (colname in names(data_df)){
  data[[colname]] <- data_df[[colname]]
}



results <- jags(
  data = data,
  # inits = initial_values,
  parameters.to.save = c("OR", "theta"),
  model.file = model_code,
  n.chains = 1, # length(initial_values),
  # n.adapt = 100,
  n.iter = 50000,
  n.burnin = 20000 #, n.thin = 2
)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 34
##    Unobserved stochastic nodes: 1
##    Total graph size: 108
## 
## Initializing model
## 
## Adaptive phase..... 
## Adaptive phase complete 
##  
## 
##  Burn-in phase, 20000 iterations x 1 chains 
##  
## 
## Sampling from joint posterior, 30000 iterations x 1 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.
summary(results)
##                mean         sd       2.5%        25%        50%        75%
## OR        0.8450758 0.02730002  0.7929252  0.8263773  0.8446826  0.8633376
## theta    -0.1688505 0.03229762 -0.2320263 -0.1907038 -0.1687943 -0.1469495
## deviance 83.6143430 1.40159071 82.6188310 82.7223690 83.0754879 83.9443157
##               97.5% Rhat n.eff overlap0 f
## OR        0.8996653   NA     1        0 1
## theta    -0.1057325   NA     1        0 1
## deviance 87.5968764   NA     1        0 1
plot(results)

Cholesterol_fixed.txt

This is the file for Model 4.1 from the Exercises, referenced on p. 89.

Model 1: Generic fixed effect model applied to the cholesterol lowering dataset

model_code <- "model 
{ 
  for (i in 1:k)
  {
    P[i] <- 1/V[i]
    logor[i] ~ dnorm(d, P[i])
  }
    d ~ dnorm(0, 1.0E-5)
    OR <- exp(d)
}" %>% textConnection

# DATA: re-used from previous model. Here the name of the second column in the data table is `V` instead of `var`. I changed that in the previous model because the name `var` seemed to be causing problems.

# data <- list(k=34)
# 
# data_df <- "logor V
# 0.4893853 0.1987757
# -0.0683811    0.054427
# -0.4876372    0.0730886
# -0.1630902    0.1108992
# 1.469676  1.186739
# 0.0057637 0.0608936
# -0.1496961    0.0799553
# -1.94591  2.351288
# -0.0470209    0.0194635
# -0.3824268    0.0552552
# -0.5508629    0.0615112
# -0.0984401    0.0510244
# 0.2939126 0.1708592
# -0.4883528    1.582492
# -0.4354064    0.0030637
# -0.2885568    0.0496847
# 0.4696847 0.0716345
# 0.2692525 0.0101539
# -1.563976 0.8819983
# -0.3517397    0.3733933
# -0.0480808    0.0298812
# -1.109251 2.687944
# 0.0129744 0.0469846
# -0.4124823    0.0383304
# -0.0233855    0.0200697
# 0.0803459 0.0082181
# -0.2820902    0.0420423
# -1.098612 2.707904
# 0.521464  2.696409
# -0.2501722    0.477032
# 1.025401  0.3643915
# -0.7528253    0.0676265
# -1.580711 2.272708
# 0.5030903 0.1426682" %>% read.delim(text=., sep="\t")
# 
# for (colname in names(data_df)){
#   data[[colname]] <- data_df[[colname]]
# }

# INITIAL VALUES: This is a list of lists, one per chain.
initial_values <- list(
  list(d=0)
)

results <- jags(
  data = data,
  inits = initial_values,
  parameters.to.save = c("OR", "d"),
  model.file = model_code,
  n.chains = length(initial_values),
  n.adapt = 100,
  n.iter = 50000,
  n.burnin = 20000,
  n.thin = 2
)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 34
##    Unobserved stochastic nodes: 1
##    Total graph size: 108
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 1 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 20000 iterations x 1 chains 
##  
## 
## Sampling from joint posterior, 30000 iterations x 1 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.
summary(results)
##                mean         sd       2.5%        25%        50%        75%
## OR        0.8447898 0.02718736  0.7929180  0.8262318  0.8446778  0.8628069
## d        -0.1691850 0.03217537 -0.2320354 -0.1908800 -0.1688000 -0.1475643
## deviance 83.6067409 1.39355862 82.6187542 82.7180516 83.0653805 83.9429405
##               97.5% Rhat n.eff overlap0 f
## OR        0.8981821   NA     1        0 1
## d        -0.1073824   NA     1        0 1
## deviance 87.3827943   NA     1        0 1
plot(results)

Cholesterol_random.txt

This is the file for Model 4.2 from the Exercises, referenced on p. 90.

Model 2: Generic random effect model applied to the cholesterol lowering dataset

model_code <- "model 
{ 
  for (i in 1:k)
    {
      P[i] <- 1/V[i]
      logor[i] ~ dnorm(delta[i], P[i])
      delta[i] ~ dnorm(d, prec) 
    }   
    d ~ dnorm(0, 1.0E-5)
    OR <- exp(d)
    
    tau ~ dunif(0,10)
    tau.sq <- tau*tau
    prec <- 1/(tau.sq)
    d.new ~ dnorm(d, prec)
}" %>% textConnection

# DATA: re-used from previous chunks.


initial_values <- list(
  list(d=0, 
       delta=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), 
       d.new=0)
)

results <- jags(
  data = data,
  inits = initial_values,
  parameters.to.save = c("OR", "d"),
  model.file = model_code,
  n.chains = length(initial_values),
  n.adapt = 100,
  n.iter = 50000,
  n.burnin = 20000,
  n.thin = 2
)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 34
##    Unobserved stochastic nodes: 37
##    Total graph size: 147
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 1 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 20000 iterations x 1 chains 
##  
## 
## Sampling from joint posterior, 30000 iterations x 1 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.
summary(results)
##                mean         sd       2.5%        25%        50%         75%
## OR        0.8912361 0.06032464  0.7768036  0.8508168  0.8892777  0.92939818
## d        -0.1174281 0.06755136 -0.2525678 -0.1615585 -0.1173457 -0.07321802
## deviance 26.8779939 6.35016609 15.4871660 22.4505696 26.5334845 30.85735636
##               97.5% Rhat n.eff overlap0         f
## OR        1.0157291   NA     1        0 1.0000000
## d         0.0156067   NA     1        1 0.9590667
## deviance 40.3147275   NA     1        0 1.0000000

Cholesterol_random_logistic.txt

This is the file for Model 4.3 in the Exercises, p90.

model_code <- "model
    {
       for( i in 1 : k ) {
        rA[i] ~ dbin(pA[i], nA[i])
           rB[i] ~ dbin(pB[i], nB[i])
        logit(pA[i]) <- mu[i]
        logit(pB[i]) <- mu[i] + delta[i]
        mu[i] ~ dnorm(0.0,1.0E-5)
        delta[i] ~ dnorm(d, prec)
       }
        OR <- exp(d)
        d ~ dnorm(0.0,1.0E-6)
    
        tau~dunif(0,10)
        tau.sq<-tau*tau
        prec<-1/(tau.sq)
    }

" %>% textConnection

# DATA
data <- list(k=34)

data_df <- "nB  nA  rB  rA
50  50  17  12
285 147 70  38
156 119 37  40
123 129 20  24
54  26  8   1
427 143 81  27
199 194 28  31
30  33  0   3
424 422 174 178
206 206 41  55
244 253 31  51
350 367 42  48
47  48  23  20
23  29  1   2
5552    2789    1025    723
1149    1129    37  48
221 237 39  28
5331    5296    236 181
88  30  2   3
71  72  5   7
1906    1900    68  71
94  94  0   1
2051    2030    44  43
279 276 61  82
1018    1015    111 113
4541    4516    269 248
421 417 49  62
48  49  0   1
94  52  1   0
79  78  4   5
6582    1663    33  3
204 202 28  51
30  60  0   4
311 317 19  12" %>% read.delim(text=., sep="\t")

for (colname in names(data_df)){
  data[[colname]] <- data_df[[colname]]
}

initial_values = list(
  list(d=0, 
       tau=1, 
       delta=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), # rep(0, data$k)
       mu=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))
  )

results <- jags(
  data = data,
  inits = initial_values,
  parameters.to.save = c("OR", "theta"),
  model.file = model_code,
  n.chains = length(initial_values),
  n.adapt = 100,
  n.iter = 50000,
  n.burnin = 20000 #, n.thin = 2
)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 68
##    Unobserved stochastic nodes: 70
##    Total graph size: 317
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 1 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 20000 iterations x 1 chains 
##  
## 
## Sampling from joint posterior, 30000 iterations x 1 chains 
## 
## Warning in FUN(X[[i]], ...): Failed to set trace monitor for theta
## Variable theta not found
## 
## Calculating statistics....... 
## 
## Done.
summary(results)
##                 mean          sd        2.5%         25%         50%
## OR         0.8931273  0.06185544   0.7784156   0.8515065   0.8905472
## deviance 385.6625551 10.82747134 365.8803706 378.0923377 385.1371885
##                  75%      97.5% Rhat n.eff overlap0 f
## OR         0.9320761   1.022215   NA     1        0 1
## deviance 392.6321371 408.210312   NA     1        0 1
plot(results)

Cholesterol_random_logistic_2

Also referenced on p90.

Model 3 (Cont) Question 5: Direct OR model applied to the cholesterol lowering dataset

model_code <- "

model
{
  for( i in 1 : k ) {
    rA[i] ~ dbin(pA[i], nA[i])
    rB[i] ~ dbin(pB[i], nB[i])
    logit(pA[i]) <- mu[i]
    logit(pB[i]) <- mu[i] + delta[i]
    mu[i] ~ dnorm(0.0,1.0E-5)
    delta[i] ~ dnorm(d, prec)
  }
  OR <- exp(d)
  d ~ dnorm(0.0,1.0E-6)
  
  tau~dunif(0,10)
  tau.sq<-tau*tau
  prec<-1/(tau.sq)
  
  # Storing theta.pooled in delta[37] to add to caterpillar plot
  delta[37] <- d
  
  # Replicate from a predictive distribution for a future study
  delta[39] ~ dnorm(d, prec)
  
}
" %>% textConnection

# DATA
data <- list(k=34)

data_df <- "nB  nA  rB  rA
50  50  17  12
285 147 70  38
156 119 37  40
123 129 20  24
54  26  8   1
427 143 81  27
199 194 28  31
30  33  0   3
424 422 174 178
206 206 41  55
244 253 31  51
350 367 42  48
47  48  23  20
23  29  1   2
5552    2789    1025    723
1149    1129    37  48
221 237 39  28
5331    5296    236 181
88  30  2   3
71  72  5   7
1906    1900    68  71
94  94  0   1
2051    2030    44  43
279 276 61  82
1018    1015    111 113
4541    4516    269 248
421 417 49  62
48  49  0   1
94  52  1   0
79  78  4   5
6582    1663    33  3
204 202 28  51
30  60  0   4
311 317 19  12" %>% read.delim(text=., sep="\t")

for (colname in names(data_df)){
  data[[colname]] <- data_df[[colname]]
}

initial_values = list(
  list(
    d=0, 
    tau=2, 
    delta=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,NA,NA,NA,NA,1),  
    mu=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))

)


results <- jags(
  data = data,
  inits = initial_values,
  parameters.to.save = c("OR", "theta"),
  model.file = model_code,
  n.chains = length(initial_values),
  n.adapt = 100,
  n.iter = 50000,
  n.burnin = 20000 #, n.thin = 2
)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 68
##    Unobserved stochastic nodes: 71
##    Total graph size: 318
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 1 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 20000 iterations x 1 chains 
##  
## 
## Sampling from joint posterior, 30000 iterations x 1 chains 
## 
## Warning in FUN(X[[i]], ...): Failed to set trace monitor for theta
## Variable theta not found
## 
## Calculating statistics....... 
## 
## Done.
summary(results)
##                 mean          sd        2.5%        25%         50%        75%
## OR         0.8944183  0.06199873   0.7784412   0.853385   0.8918417   0.933044
## deviance 385.5913474 10.71695257 366.1612898 378.094530 385.0791317 392.534913
##               97.5% Rhat n.eff overlap0 f
## OR         1.026008   NA     1        0 1
## deviance 408.177713   NA     1        0 1
plot(results)

Cholesterol_prior_sens

This file is referenced for Model 4.4, p91.

Model 4: Direct OR model with built in sensitivity analysis to prior

RUNTIME ERROR:
Compilation error on line 14.
Attempt to redefine node rBx[1,1]

Original WinBUGS code:

model
{
    
  ## replicating the data 2 times to fit to 3 models in total where j indexes alternative prior distributions

  for (j in 1:3)
  {
     
    for( i in 1 : k )
    {
      rAx[j,i] ~ dbin(pA[j,i], nA[i])
      rAx[j,i] <- rA[i]
      rBx[j,i] ~ dbin(pB[j,i], nB[i])
      rBx[j,i] <- rB[i]
      logit(pA[j,i]) <- mu[j,i]
      logit(pB[j,i]) <- mu[j,i] + delta[j,i]
      mu[j,i] ~ dnorm(0.0,1.0E-5)
      delta[j,i] ~ dnorm(d[j], prec[j])
    }
  }

  for (j in 1:3)
  {
    OR[j] <- exp(d[j])
    d[j] ~ dnorm(0.0,1.0E-6)
    m[j] ~ dnorm(0.0,1.0E-6)
  }

    
  # Setting three different priors for the between study variance
  
  prec[1] <-1/tau.sq[1]
  tau.sq[1] <-tau[1]*tau[1]
  tau[1] ~ dnorm(0,1.0E-6)I(0,)
  
  prec[2] ~ dgamma(0.001,0.001)
  tau.sq[2] <- 1/prec[2]
  tau[2] <- sqrt(tau.sq[2])
  
  tau[3] ~ dunif(0,50)
  tau.sq[3] <- tau[3]*tau[3]
  prec[3] <- 1/tau.sq[3]
    
}
model_code <- "
data{
  for (j in 1:3)
  {
    for( i in 1 : k )
    {
      rAx[j,i] <- rA[i]
      rBx[j,i] <- rB[i]
    }
  }

}

model
{
    
  ## replicating the data 2 times to fit to 3 models in total where j indexes alternative prior distributions

  for (j in 1:3)
  {
     
    for( i in 1 : k )
    {
      rAx[j,i] ~ dbin(pA[j,i], nA[i])
      # rAx[j,i] <- rA[i]   # !!! Attempt to redefine node
      rBx[j,i] ~ dbin(pB[j,i], nB[i])
      # rBx[j,i] <- rB[i]   # !!! Attempt to redefine node
      logit(pA[j,i]) <- mu[j,i]
      logit(pB[j,i]) <- mu[j,i] + delta[j,i]
      mu[j,i] ~ dnorm(0.0,1.0E-5)
      delta[j,i] ~ dnorm(d[j], prec[j])
    }
  }

  for (j in 1:3)
  {
    OR[j] <- exp(d[j])
    d[j] ~ dnorm(0.0,1.0E-6)
    m[j] ~ dnorm(0.0,1.0E-6)
  }

    
  # Setting three different priors for the between study variance
  
  prec[1] <-1/tau.sq[1]
  tau.sq[1] <-tau[1]*tau[1]
  tau[1] ~ dnorm(0,1.0E-6)I(0,)
  
  prec[2] ~ dgamma(0.001,0.001)
  tau.sq[2] <- 1/prec[2]
  tau[2] <- sqrt(tau.sq[2])
  
  tau[3] ~ dunif(0,50)
  tau.sq[3] <- tau[3]*tau[3]
  prec[3] <- 1/tau.sq[3]
    
}" %>% textConnection

# DATA
data <- list(k=34)

data_df <- "nB  nA  rB  rA
50  50  17  12
285 147 70  38
156 119 37  40
123 129 20  24
54  26  8   1
427 143 81  27
199 194 28  31
30  33  0   3
424 422 174 178
206 206 41  55
244 253 31  51
350 367 42  48
47  48  23  20
23  29  1   2
5552    2789    1025    723
1149    1129    37  48
221 237 39  28
5331    5296    236 181
88  30  2   3
71  72  5   7
1906    1900    68  71
94  94  0   1
2051    2030    44  43
279 276 61  82
1018    1015    111 113
4541    4516    269 248
421 417 49  62
48  49  0   1
94  52  1   0
79  78  4   5
6582    1663    33  3
204 202 28  51
30  60  0   4
311 317 19  12" %>% read.delim(text=., sep="\t")

for (colname in names(data_df)){
  data[[colname]] <- data_df[[colname]]
}

initial_values = list(

  list(
    prec=c(NA,0.1,NA),  
    tau=c(0.5, NA, 0.5), 
    mu = structure(.Data=c(0,0,0,0,0,0,0,0,0,0,  0,0,0,0,0,0,0,0,0,0,  0,0,0,0,0,0,0,0,0,0,   0,0,0,0,
0,0,0,0,0,0,0,0,0,0,  0,0,0,0,0,0,0,0,0,0,  0,0,0,0,0,0,0,0,0,0,   0,0,0,0,
0,0,0,0,0,0,0,0,0,0,  0,0,0,0,0,0,0,0,0,0,  0,0,0,0,0,0,0,0,0,0,   0,0,0,0),.Dim=c(3,34)),

    delta = structure(.Data=c(
0,0,0,0,0,0,0,0,0,0,  0,0,0,0,0,0,0,0,0,0,  0,0,0,0,0,0,0,0,0,0,   0,0,0,0,
0,0,0,0,0,0,0,0,0,0,  0,0,0,0,0,0,0,0,0,0,  0,0,0,0,0,0,0,0,0,0,   0,0,0,0,
0,0,0,0,0,0,0,0,0,0,  0,0,0,0,0,0,0,0,0,0,  0,0,0,0,0,0,0,0,0,0,   0,0,0,0),.Dim=c(3,34)) #)
    , 
    d =c(0.1,0.1,0.1) 
  )

)


results <- jags(
  data = data,
  inits = initial_values,
  parameters.to.save = c("OR", "d", "theta"),
  model.file = model_code,
  n.chains = length(initial_values),
  n.adapt = 100,
  n.iter = 50000,
  n.burnin = 20000 #, n.thin = 2
)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling data graph
##    Resolving undeclared variables
##    Allocating nodes
##    Initializing
##    Reading data back into data table
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 204
##    Unobserved stochastic nodes: 213
##    Total graph size: 875
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 1 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 20000 iterations x 1 chains 
##  
## 
## Sampling from joint posterior, 30000 iterations x 1 chains 
## 
## Warning in FUN(X[[i]], ...): Failed to set trace monitor for theta
## Variable theta not found
## 
## Calculating statistics....... 
## 
## Done.
summary(results)
##                  mean          sd         2.5%          25%          50%
## OR[1]       0.8942644  0.06234855    0.7779298    0.8527763    0.8915863
## OR[2]       0.8927113  0.05906542    0.7816431    0.8530278    0.8910338
## OR[3]       0.8930435  0.06182829    0.7770039    0.8514910    0.8910776
## d[1]       -0.1141740  0.06955541   -0.2511190   -0.1592580   -0.1147531
## d[2]       -0.1156736  0.06604820   -0.2463570   -0.1589631   -0.1153729
## d[3]       -0.1155074  0.06908539   -0.2523098   -0.1607664   -0.1153238
## deviance 1157.6300866 18.61739198 1122.9159571 1144.6421859 1157.0548490
##                    75%        97.5% Rhat n.eff overlap0         f
## OR[1]       0.93307483 1.024227e+00   NA     1        0 1.0000000
## OR[2]       0.92989391 1.015263e+00   NA     1        0 1.0000000
## OR[3]       0.93218913 1.019841e+00   NA     1        0 1.0000000
## d[1]       -0.06926988 2.393851e-02   NA     1        1 0.9502667
## d[2]       -0.07268477 1.514736e-02   NA     1        1 0.9590667
## d[3]       -0.07021956 1.964635e-02   NA     1        1 0.9543667
## deviance 1169.98851875 1.195988e+03   NA     1        0 1.0000000
plot(results)